library(ggplot2)
library(reshape2)
library(tidyverse)
library(dplyr)
library(e1071) # asimetria
library(plyr)
library(fitdistrplus) # grafica Cullen-Frey
library(dslabs) # datos PCA
library(factoextra) # graficas PCA
library(purrr)
library(bestNormalize)
library(car)
library(ggiraph)
library(predict3d)
library(grid)
library(gridExtra)
require(nortest)

Limpieza de datos

Abrimos el archivo y observamos los campos que se encuentran en nuestros datos.

# Abrir el archivo de datos
WCGS <- read.csv("source/wcgs.csv")
head(WCGS)
##   age arcus behpat      bmi chd69 chol dbp dibpat height    id    lnsbp
## 1  50     1     A1 31.32101    No  249  90 Type A     67  2343 4.882802
## 2  51     0     A1 25.32858    No  194  74 Type A     73  3656 4.787492
## 3  59     1     A1 28.69388    No  258  94 Type A     70  3526 5.062595
## 4  51     1     A1 22.14871    No  173  80 Type A     69 22057 4.836282
## 5  44     0     A1 22.31303    No  214  80 Type A     71 12927 4.836282
## 6  47     0     A1 27.11768    No  206  76 Type A     64 16029 4.753590
##     lnwght ncigs sbp smoke         t1 time169 typchd69       uni weight wghtcat
## 1 5.298317    25 132   Yes -1.6333529    1367        0 0.4860738    200 170-200
## 2 5.257495    25 120   Yes -4.0633659    2991        0 0.1859543    192 170-200
## 3 5.298317     0 158    No  0.6397287    2960        0 0.7277991    200 170-200
## 4 5.010635     0 126    No  1.1217681    3069        0 0.6244636    150 140-170
## 5 5.075174     0 126    No  2.4250107    3081        0 0.3789776    160 140-170
## 6 5.062595    80 116   Yes -0.7875199    2114        0 0.7355005    158 140-170
##    agec
## 1 46-50
## 2 51-55
## 3 56-60
## 4 51-55
## 5 41-45
## 6 46-50

En contramos que los campos son:

  • id Subject ID:

  • age0 Age: age in years

  • height0 Height: height in inches

  • weight0 Weight: weight in pounds

  • sbp0 Systolic blood pressure: mm Hg

  • dbp0 Diastolic blood pressure: mm Hg

  • chol0 Cholesterol: mg/100 ml

  • behpat0 Behavior pattern:

  • ncigs0 Smoking: Cigarettes/day

  • dibpat0 Dichotomous behavior pattern: 0 = Type B; 1 = Type A

  • chd69 Coronary heart disease event: 0 = none; 1 = yes

  • typechd to be done

  • time169 Observation (follow up) time: Days

  • arcus0 Corneal arcus: 0 = none; 1 = yes

Datos faltantes

Ante estos datos primero que nada deberemos realizar una limpieza de los datos faltantes, esto con el propósito de evitar que se genere ruido en nuestros datos.

Para esto haremos un comparativo entre la longitud del dataframe antes de omitir datos faltantes y despues.

# Datos antes de omitir faltantes
length(WCGS$age)
## [1] 3154
# Datos despues de eliminar faltantes
WCGS <- na.omit(WCGS)
length(WCGS$age)
## [1] 3101

Una vez que realizamos el comparativo podemos observar que se perdieron un total de 53 datos que faltaban en el dataframe.

Limpieza de valores extremos

En este caso realizaremos un boxplot de las multiples variables cuantitativas en los datos, estas variables seran:

  • Colesterol
  • Peso
  • Edad
  • Altura
  • Presión sístolica
  • Presión diastólica
# Matriz para el layout
x <- c(1, 2, 3) # tres gráficos diferentes
m <- matrix(x, ncol = 2)

nf <- layout(m)

# Plot de variable respuesta
Plot1 <- ggplot(WCGS, aes(x = chd69, y = weight)) + 
  geom_boxplot(colour=c("green","blue")) +
  labs(title = "Boxplot del peso relacionado con la \nenfermedad coronaria ", x = "Enfermedad Coronaria", y = "Peso") + theme(plot.title = element_text(size=7.5))

# Plot de variable predictora 1
Plot2 <- ggplot(WCGS, aes(x = chd69, y = height)) + 
  geom_boxplot(colour=c("green","blue")) +
  labs(title = "Boxplot de la altura relacionada con la \nenfermedad coronaria ", x = "Enfermedad Coronaria", y = "Altura") + theme(plot.title = element_text(size=7.5))

# Plot de variable predictora 2
Plot3 <- ggplot(WCGS, aes(x = chd69, y = chol)) + 
  geom_boxplot(colour=c("green","blue")) +
  labs(title = "Boxplot del colesterol relacionado con la \nenfermedad coronaria ", x = "Enfermedad Coronaria", y = "Colesterol") 



Plots <- list(Plot1,Plot2,Plot3)

layout <- rbind(c(1,2),c(3))

grid.arrange(grobs = Plots,layout_matrix = layout)

De acuerdo con la información del estudio las medidas se realizaron mediante el uso del sistema ingles, por este motivo tendremos que darle una interpretación que nos arrojan estos datos.

  • Peso: En este caso observamos que las personas con mayor peso se encuentran en una medición de 320 libras y los de menor en 100 libras, lo cual es equivalente a 145kg y 45kg respectivamente. De modo que, en este campo consideramos que si bien hay personas con bajo peso y personas obesas, consideramos que los datos no son tan extremos como para descartarlo.

  • Altura: En este caso observamos que las personas con mayor altura se encuentran en una medición de 79 pulgadas y los de menor en 60 pulgadas, lo cual es equivalente a 2 metros y 1.52 metros. De modo que, en este campo consideramos que la altura no presenta ruido en los datos.

  • Colesterol: En este caso encontramos personas en las que se anoto una cantidad de colesterol muy grande cercana a 600 por este motivo consideramos que lo mas adecuado es retirar estos datos, estableciendo un punto de corte en 400 que es donde observamos el agrupamiento.

# Matriz para el layout
x <- c(1, 2, 3) # tres gráficos diferentes
m <- matrix(x, ncol = 2)

nf <- layout(m)

# Plot de variable predictora 3
Plot4 <- ggplot(WCGS, aes(x = chd69, y = sbp)) + 
  geom_boxplot(colour=c("green","blue")) +
  labs(title = "Boxplot de la P.sistólica relacionada con la \nenfermedad coronaria ", x = "Enfermedad Coronaria", y = "Presión sistolica") + theme(plot.title = element_text(size=7.5))

# Plot de variable predictora 3
Plot5 <- ggplot(WCGS, aes(x = chd69, y = dbp)) + 
  geom_boxplot(colour=c("green","blue")) +
  labs(title = "Boxplot de la P.diastólica relacionada \ncon enfermedad coronaria ", x = "Enfermedad Coronaria", y = "Presión diastolica") + theme(plot.title = element_text(size=7.5))

# Plot de variable predictora 3
Plot6 <- ggplot(WCGS, aes(x = chd69, y = age)) + 
  geom_boxplot(colour=c("green","blue")) +
  labs(title = "Boxplot de la edad relacionada \ncon enfermedad coronaria ", x = "Enfermedad Coronaria", y = "Edad") 

Plots <- list(Plot4,Plot5,Plot6)

layout <- rbind(c(1,2),c(3))

grid.arrange(grobs = Plots,layout_matrix = layout)

En este apartado podemos observar que lo indicado por el articulo es concordante con los datos pues todos los pacientes anotados se encuentran en un rango de edad de 39 a 60 años. Por otra parte se encontraron valores altos y bajos en las presiones sistolica y diastolica, pero no fueron tan extremos como para considerar limpiarlos.

Eliminamos los datos extremos de colesterol.

newChol <- (WCGS$chol <= 400)
WCGS <- WCGS[newChol,]
length(WCGS$age)
## [1] 3099

Visualización de datos

Values <- data.frame(Edad = as.numeric(WCGS$age), IMC = as.numeric(WCGS$bmi),Altura = as.numeric(WCGS$height), Coronaria = WCGS$chd69, Colesterol = as.numeric(WCGS$chol))
head(Values)
##   Edad      IMC Altura Coronaria Colesterol
## 1   50 31.32101     67        No        249
## 2   51 25.32858     73        No        194
## 3   59 28.69388     70        No        258
## 4   51 22.14871     69        No        173
## 5   44 22.31303     71        No        214
## 6   47 27.11768     64        No        206
max(Values$IMC)
## [1] 38.94737
min(Values$IMC)
## [1] 11.19061

Categorizar datos

Para unas graficas mas entendibles categorizaremos los datos colocados en nuestro nuevo data frame.

################## Categorizar pesos #################
obex <- Values["IMC"] > 39.9
obeso <- (Values["IMC"] > 29.9) & (Values["IMC"] <= 39.9)
sobrep <- (Values["IMC"] > 24.9) & (Values["IMC"] <= 29.9)
Norm <- (Values["IMC"] > 18.5) & (Values["IMC"] <= 24.9)
bajo <- Values["IMC"] <= 18.5

Values["IMC"][obex] <- "Extremo"
Values["IMC"][obeso] <- "obeso"
Values["IMC"][sobrep] <- "Sobrepeso"
Values["IMC"][Norm] <- "Normal"
Values["IMC"][bajo] <- "Bajo"

################## Categorizar altura #################
Bajo <- Values["Altura"] <= 65
Promedio <- (Values["Altura"] > 65) & (Values["Altura"] < 71)
Alto <- Values["Altura"] >= 71

Values["Altura"][Bajo] <- "Bajo"
Values["Altura"][Promedio] <- "Normal"
Values["Altura"][Alto] <- "Alto"

################## Categorizar edad #################
Adulto <- Values["Edad"] <= 45
Mayor <- (Values["Edad"] > 45) & (Values["Edad"] < 55)
Viejo <- Values["Edad"] >= 55

Values["Edad"][Adulto] <- "Adulto"
Values["Edad"][Mayor] <- "Adulto mayor"
Values["Edad"][Viejo] <- "Viejo"

################ Categorizar colesterol ###################
CNorm <- Values["Colesterol"] <= 200
CNAlt <- (Values["Colesterol"] > 200) & (Values["Colesterol"] < 240)
CAlt <- Values["Colesterol"] >= 240

Values["Colesterol"][CNorm] <- "Normal"
Values["Colesterol"][CNAlt] <- "Normal-Alto"
Values["Colesterol"][CAlt] <- "Alto"

head(Values)
##           Edad       IMC Altura Coronaria  Colesterol
## 1 Adulto mayor     obeso Normal        No        Alto
## 2 Adulto mayor Sobrepeso   Alto        No      Normal
## 3        Viejo Sobrepeso Normal        No        Alto
## 4 Adulto mayor    Normal Normal        No      Normal
## 5       Adulto    Normal   Alto        No Normal-Alto
## 6 Adulto mayor Sobrepeso   Bajo        No Normal-Alto
# Realizar conteo de categorias
table(Values$IMC)
## 
##      Bajo    Normal     obeso Sobrepeso 
##        21      1770        86      1222
IMC <- table(Values$IMC)
IMC <- as.data.frame(IMC)

table(Values$Colesterol)
## 
##        Alto      Normal Normal-Alto 
##        1087         842        1170
Colesterol <- table(Values$Colesterol)
Colesterol <- as.data.frame(Colesterol)

table(Values$Altura)
## 
##   Alto   Bajo Normal 
##   1204    132   1763
Altura <- table(Values$Altura)
Altura <- as.data.frame(Altura)

table(Values$Edad)
## 
##       Adulto Adulto mayor        Viejo 
##         1605         1182          312
Edad <- table(Values$Edad)
Edad <- as.data.frame(Edad)

Graficar categorias

Por medio de recursos graficos observaremos como se encuentran distribuidos los pacientes del dataset de acuerdo a sus caracteristicas.

# Matriz para el layout
x <- c(1, 2) # cuatro gráficos diferentes
m <- matrix(x, ncol = 2)

nf <- layout(m)

# Plot de variable respuesta
Plot1 <- ggplot(IMC, aes(x = Var1, y = Freq, fill = c("Bajo","Normal", "Obeso","Sobrepeso"))) +
  geom_bar(stat = "identity") + labs(title = "Número de Pacientes pertenecientes a cada \ngrupo de IMC", x = "Nivel de IMC", y = "Número de pacientes") + theme(axis.ticks.x = element_blank(), axis.text.x=element_blank()) + guides(fill=guide_legend(title="Nivel IMC")) + theme(legend.key.size = unit(0.2,"line"))

# Plot de variable predictora 1
Plot2 <- ggplot(Colesterol, aes(x = Var1, y = Freq, fill = c("Alto","Normal", "Normal-Alto"))) +
  geom_bar(stat = "identity") + labs(title = "Número de Pacientes pertenecientes a cada \ngrupo de Colesterol", x = "Nivel de Colesterol", y = "Número de pacientes") + theme(axis.ticks.x = element_blank(), axis.text.x=element_blank()) + guides(fill=guide_legend(title="Nivel Colesterol")) +  theme(legend.key.size = unit(0.2,"line"))


Plots <- list(Plot1,Plot2)

layout <- rbind(c(1,2))

grid.arrange(grobs = Plots,layout_matrix = layout)

# Matriz para el layout
x <- c(1, 2) # cuatro gráficos diferentes
m <- matrix(x, ncol = 2)

nf <- layout(m)

# Plot de variable predictora 2
Plot3 <- ggplot(Altura, aes(x = Var1, y = Freq, fill = c("Alto","Bajo","Normal"))) +
  geom_bar(stat = "identity") + labs(title = "Número de Pacientes pertenecientes a cada \ngrupo de Altura", x = "Nivel de Altura", y = "Número de pacientes") + theme(axis.ticks.x = element_blank(), axis.text.x=element_blank()) + guides(fill=guide_legend(title="Nivel Altura")) +  theme(legend.key.size = unit(0.2,"line"))

# Plot de variable predictora 3
Plot4 <- ggplot(Edad, aes(x = Var1, y = Freq, fill = c("Adulto","Adulto M", "Viejo"))) +
  geom_bar(stat = "identity") + labs(title = "Número de Pacientes pertenecientes a cada \ngrupo de Edad", x = "Nivel de Edad", y = "Número de pacientes") + theme(axis.ticks.x = element_blank(), axis.text.x=element_blank()) + guides(fill=guide_legend(title="Grupo de Edad")) +  theme(legend.key.size = unit(0.2,"line"))


Plots <- list(Plot3,Plot4)

layout <- rbind(c(1,2))

grid.arrange(grobs = Plots,layout_matrix = layout)

Ahora pensamos que el colesterol alto puede estar muy relacionado con el desarrollo de una enfermedad coronaria.

# Realizar el conteo cardiovascular dependiendo del nivel de colesterol
Group <- tapply(Values$Colesterol, 
                  list(Values$Colesterol, Values$Coronaria),
                  table)

Group <- as.data.frame(Group)
Group$Colesterol <- rownames(Group)
Group
##               No Yes  Colesterol
## Alto         951 136        Alto
## Normal       808  34      Normal
## Normal-Alto 1091  79 Normal-Alto
# Generar el formato long del dataframe
Group <- reshape(Group, #Objeto a transformar
        direction = "long", #Necesitamos formato long
        varying = list(names(Group)[1:2]), #Las columnas a alargar (las que tienen las cuentas) son la 1 y la 2
        v.names = "Counts", #Las cuentas se almacenaran en la columna "Counts"
        idvar = "Colesterol", #La variable con los identificadores del grupo de actividad fisica
        timevar = "Coronaria", #Variable que vamos a alargar, que corresponde a las columnas "no" y "yes"
        times = c("no", "yes"))
row.names(Group) <- NULL
Group
##    Colesterol Coronaria Counts
## 1        Alto        no    951
## 2      Normal        no    808
## 3 Normal-Alto        no   1091
## 4        Alto       yes    136
## 5      Normal       yes     34
## 6 Normal-Alto       yes     79
# Barplot que refleja los niveles de colesterol dependiendo el grupo
ggplot(Group, aes(x = Colesterol, y = Counts, fill = Coronaria)) +
  geom_bar(stat = "identity", position=position_dodge()) + labs(title = "Número de Pacientes pertenecientes a cada \ngrupo de colesterol y enfermedad cardiovascular", x = "Nivel de colesterol", y = "Número de pacientes") + scale_fill_discrete(name = "Cardiovascular", labels = c("Ausencia", "Presencia"))

Análisis de la varianza

Debido a lo anterior, hemos decidido analizar si la diferencia del nivel de colesterol, entre pacientes sanos y enfermos, es significativa.

# Importar la libreria requerida para realizar la codificacion de datos
library(qdap)
# Establecer un tipo de variable valido para la categoria de genero, para no tener errores
WCGS$chd69 <- multigsub(sort(unique(WCGS$chd69)), c("Ausencia", "Presencia"), WCGS$chd69)
table(WCGS$chd69)
## 
##  Ausencia Presencia 
##      2850       249

Como se puede observar, los dos conjuntos (sanos y enfermos) tienen tamaños distintos. A continuación, analizaremos la distribución de cada uno.

# Calcular el promedio de cada conjunto
aggregate(chol~chd69, data = WCGS, FUN = mean)
##       chd69     chol
## 1  Ausencia 224.1846
## 2 Presencia 247.9880

En promedio, los pacientes sanos tienen menor nivel de colesterol que los pacientes enfermos, en este conjunto de datos.

# Calcular la desviacion estandar de cada conjunto
aggregate(chol~chd69, data = WCGS, FUN = sd)
##       chd69     chol
## 1  Ausencia 42.25111
## 2 Presencia 41.89460

La desviación estándar en ambos conjuntos es similar, lo cual podría darnos indicios de que comparten la misma distribución. A continuación, aplicaremos el test Lilliefors para evaluar de forma analítica la normalidad en ambos conjuntos. Escogimos esta prueba ya que es la alternativa a Shapiro-Wilk cuando el número de observaciones es mayor de 50, como es el caso. La Hipótesis nula (H0) de esta prueba es que la muestra tiene una distribución normal, y la hipótesis alternativa (HA) sería que las muestras no tienen una distribución normal. Por lo tanto, si el p-valor es menor a 0.05 se puede rechazar la H0 y aceptar la hipótesis alternativa.

# Separar los datos por genero en estructuras distintas
healthy <- WCGS[WCGS$chd69=="Ausencia",]
sick <- WCGS[WCGS$chd69=="Presencia",]
# Comprobar la normalidad de nuestros datos
lillie.test(healthy$chol)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  healthy$chol
## D = 0.037532, p-value = 6.311e-10
lillie.test(sick$chol)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  sick$chol
## D = 0.048163, p-value = 0.1715

El p-valor de la prueba Lillefors aplicada en el conjunto de individuos sanos es menor a 0.05, lo cual indica que su distribución definitivamente no es normal. Sin embargo, el p-valor de la misma prueba, pero aplicada en el conjunto de pacientes enfermos, es mayor a 0.05, y por ello podemos asegurar que los datos siguen una distribución normal. Podemos dibujar las gráficas cuantil-cualtil de ambos conjuntos para comprobar lo anterior.

# Matriz para el layout
x <- c(1, 2) # dos gráficos diferentes
m <- matrix(x, ncol = 2)
nf <- layout(m)
Plot1 <- ggplot(healthy, aes(sample = chol)) + 
  stat_qq(alpha = 0.5) + 
  stat_qq_line(color = "red") + 
  labs(title = "Colesterol en pacientes sanos", x = "Teoréticos", y = "Muestra")
Plot2 <- ggplot(sick, aes(sample = chol)) + 
  stat_qq(alpha = 0.5) + 
  stat_qq_line(color = "blue") + 
  labs(title = "Colesterol en pacientes enfermos", x = "Teoréticos", y = "Muestra")
Plots <- list(Plot1,Plot2)
layout <- rbind(c(1,2))
grid.arrange(grobs = Plots,layout_matrix = layout)

Como podemos observar, el nivel de colesterol en los pacientes sanos no siguen una distribución normal, y por este motivo debemos normalizarlos.

# Normalizar los datos
set.seed(42)
healthy$NormChol <- bestNormalize(healthy$chol)$x.t
# Comprobar la normalidad de nuestros datos
lillie.test(healthy$NormChol)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  healthy$NormChol
## D = 0.018767, p-value = 0.02169

De nuevo los valores del colesterol en pacientes sanos no siguen una distribución normal, pues el p-valor de la prueba Lilliefors aplicada en los valores normalizados sigue siendo menor a 0.05. Podemos comparar el cambio entre los valores originales y normalizados en una gráfica cuantil-cuantil.

# Matriz para el layout
x <- c(1, 2) # dos gráficos diferentes
m <- matrix(x, ncol = 2)
nf <- layout(m)
Plot1 <- ggplot(healthy, aes(sample = chol)) + 
  stat_qq(alpha = 0.5) + 
  stat_qq_line(color = "red") + 
  labs(title = "Colesterol original", x = "Teoréticos", y = "Muestra")
Plot2 <- ggplot(healthy, aes(sample = NormChol)) + 
  stat_qq(alpha = 0.5) + 
  stat_qq_line(color = "blue") + 
  labs(title = "Colesterol normalizado", x = "Teoréticos", y = "Muestra")
Plots <- list(Plot1,Plot2)
layout <- rbind(c(1,2))
grid.arrange(grobs = Plots,layout_matrix = layout)

Los resultados anteriores nos indican que si bien hubo una mejora al normalizar los datos, no fue suficiente para que estos alcanzaran una distribución normal. Este hecho influye profundamente en la decisión de cuál test debe usarse para analizar la distribución de la varianza, ya que como uno de los conjuntos que se comparan no es de tipo normal, es recomendable recurrir a un test que compare la mediana de la varianza. Por lo tanto, decidí usar el test de Fligner-Killeen, ya que es un test no paramétrico que compara las varianzas basándose en la mediana y es una alternativa cuando no se cumple la condición de normalidad en las muestras.

# Comprobar la varianza constante entre grupos (homocedasticidad)
fligner.test(chol~chd69, WCGS)
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  chol by chd69
## Fligner-Killeen:med chi-squared = 0.083435, df = 1, p-value = 0.7727

El test de Fligner-Killeen considera como hipótesis nula que la varianza es igual entre los grupos y como hipótesis alternativa que no lo es. El p-value que obtuve es mayor que 0.05, por lo que la varianza entre los conjuntos es igual y por lo tanto se cumple la homocedasticidad. Para analizar la diferencia de la varianza entre ambos conjuntos, decidí usar el test Wilcoxon rank-sum, el cual es un test no paramétrico que contrasta si dos muestras proceden de poblaciones equidistribuidas. Elegí esta prueba ya que los datos cumplen con todas sus condiciones:

- Los datos son independientes.
- Los datos del nivel de colesterol son ordinales.
- El test acepta muestras que no se distribuyen de forma normal.
- Se cumple con la homocedasticidad.

Se empleará la función wilcox_test() del paquete coin, ya que los tamaños muestrales son mayores de 20.

# Importar la libreria requerida para realizar el analisis de la varianza
library(coin)
# Es necesario crear un dataframe para realizar el test
wilcox_dt <- data.frame(chol = WCGS$chol, chd69 = WCGS$chd69)
# Codificar los valores del genero en el nuevo dataframe
wilcox_dt$chd69 <- gsub("Ausencia", "N", wilcox_dt$chd69)
wilcox_dt$chd69 <- gsub("Presencia", "S", wilcox_dt$chd69)
# Aplicar el test Wilcoxon rank-sum en los datos
wilcox_test(chol ~ as.factor(chd69), data = wilcox_dt, conf.int=0.95)
## 
##  Asymptotic Wilcoxon-Mann-Whitney Test
## 
## data:  chol by as.factor(chd69) (N, S)
## Z = -8.3888, p-value < 2.2e-16
## alternative hypothesis: true mu is not equal to 0
## 95 percent confidence interval:
##  -30 -19
## sample estimates:
## difference in location 
##                    -24

Como podemos observar el p-valor de esta prueba es menor a 0.05, por lo que si existe una diferencia significativa en el nivel de colesterol entre pacientes sanos y enfermos.

Inferencia Bayesiana

Ahora analizaremos los mismos datos desde una perspectiva bayesiana.

require(BayesFactor)
## Loading required package: BayesFactor
## Loading required package: coda
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following object is masked from 'package:qdap':
## 
##     %&%
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## ************
## Welcome to BayesFactor 0.9.12-4.2. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
## 
## Type BFManual() to open the manual.
## ************
ttestBF(x = WCGS$chol[WCGS$chd69 == "Ausencia"], y = WCGS$chol[WCGS$chd69 == "Presencia"])
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 1.93261e+14 ±0%
## 
## Against denominator:
##   Null, mu1-mu2 = 0 
## ---
## Bayes factor type: BFindepSample, JZS
samples = ttestBF(x = WCGS$chol[WCGS$chd69 == "Ausencia"], y = WCGS$chol[WCGS$chd69 == "Presencia"],posterior = TRUE, iterations = 1000)

El factor calculado es 1.93261e+14 que es muy superior a 1, indicando que nuestras categorías (chol y chd69) estan altamente relacionadas. Es decir, la presencia de enfermedades coronarias tiene una gran influencia en el aumento y/o disminución del colesterol.

Regresion Lineal

En el articulo donde se realizo el estudio, se indicaba una asociacion entre la enfermedad coronaria con la conducta que presentaban los pacientes.

Normalidad IMC

Trataremos de predecir el peso de un paciente unicamente utilizando el valor del IMC con el proposito de determinar que tan preciso es el analisis de regresion con una sola variable, para esto debemos observar si hay normalidad en los datos. Para conocer esto utilizaremos una grafica cuantil-cuantil y la prueba de Shapiro.

shapiro.test(WCGS$bmi)
## 
##  Shapiro-Wilk normality test
## 
## data:  WCGS$bmi
## W = 0.9772, p-value < 2.2e-16
shapiro.test(WCGS$weight)
## 
##  Shapiro-Wilk normality test
## 
## data:  WCGS$weight
## W = 0.97801, p-value < 2.2e-16
# Matriz para el layout
x <- c(1, 2) # cuatro gráficos diferentes
m <- matrix(x, ncol = 2)

nf <- layout(m)

Plot1 <- ggplot(WCGS, aes(sample = bmi)) + 
  stat_qq(alpha = 0.5) + 
  stat_qq_line(color = "red") + 
  labs(title = "Grafica cuantil-cuantil IMC", x = "Teoréticos", y = "Muestra")

Plot2 <- ggplot(WCGS, aes(sample = weight)) + 
  stat_qq(alpha = 0.5) + 
  stat_qq_line(color = "blue") + 
  labs(title = "Grafica cuantil-cuantil Peso", x = "Teoréticos", y = "Muestra")

Plots <- list(Plot1,Plot2)

layout <- rbind(c(1,2))

grid.arrange(grobs = Plots,layout_matrix = layout)

Como podemos observar tanto el test grafico cuantil-cuantil nos indican que los datos no estan normalizados, por este motivo debemos normalizarlos. Pero, como es un modelo predictivo solo sera necesario normalizar la variable resultado.

# Normalizar los datos
set.seed(42)
bestNormalize(WCGS$weight)
## Best Normalizing transformation with 3099 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 3.5939
##  - Box-Cox: 3.5646
##  - Center+scale: 3.7046
##  - Exp(x): 83.4302
##  - Log_b(x+a): 3.5939
##  - orderNorm (ORQ): 2.9074
##  - sqrt(x + a): 3.5011
##  - Yeo-Johnson: 3.5646
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## orderNorm Transformation with 3099 nonmissing obs and ties
##  - 124 unique values 
##  - Original quantiles:
##   0%  25%  50%  75% 100% 
##   78  155  170  182  320

Una vez que ejecutamos el algoritmo podemos ver que el metodo recomendado es OrderNorm por lo que normalizamos el peso y ahora ya podemos ejecutar nuestro modelo de regresión para predecir esta variable a partir del IMC.

# Guardar datos normalizados
bestnorm <- bestNormalize(WCGS$weight)
WCGS$NormPeso <- bestnorm$x.t
head(bestnorm$x.t)
## [1]  1.4348174  1.1170894  1.4348174 -0.9784583 -0.4597549 -0.5931487
ggplot(WCGS, aes(sample = NormPeso)) + 
  stat_qq(alpha = 0.5) + 
  stat_qq_line(color = "green") + 
  labs(title = "Grafica cuantil-cuantil Peso", x = "Teoréticos", y = "Muestra")

modelo <- lm(data = WCGS, NormPeso ~ bmi)
summary(modelo)
## 
## Call:
## lm(formula = NormPeso ~ bmi, data = WCGS)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.54063 -0.39360  0.04911  0.41063  1.82148 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -7.583944   0.104571  -72.53   <2e-16 ***
## bmi          0.309329   0.004242   72.92   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.606 on 3097 degrees of freedom
## Multiple R-squared:  0.6319, Adjusted R-squared:  0.6318 
## F-statistic:  5318 on 1 and 3097 DF,  p-value: < 2.2e-16

Gracias a esto podemos establecer el modelo lineal y obtener los coeficientes de la ecuación. En este caso vemos que el valor P es mucho menor a 0.05, por lo que podemos determinar que el modelo es confiable. Los coeficientes indican que la regresión seria.

Cabe destacar que tambien se obtuvo una buena R-ajustada de 0.6318 el cual es mas cercano a 1 que a 0 y nos indica que el modelo es bueno.

\[ y = 0.309329x -7.583944 \]

Graficamos para observar el modelo de regresión con respecto a los datos.

ggplot(WCGS, aes(bmi, NormPeso)) + 
  geom_point() +
  geom_smooth(method = lm) +
  labs(title = "Grafica de dispersion entre IMC y peso con datos normalizados", x = "IMC", y = "Peso")
## `geom_smooth()` using formula 'y ~ x'

Ahora graficamos los residuales para ver que tanto se dispersan los datos de lo predicho.

WCGS$Lresiduals <- modelo$residuals
WCGS$Lpredicted <- modelo$fitted.values

ggplot(WCGS, aes(x = bmi, y = NormPeso)) +
  geom_smooth(method = lm) +  # añadimos la linea de regresion 
  geom_segment(aes(xend = bmi, yend = Lpredicted)) +  # esto añade las lineas de distancia 
  geom_point(aes(y = Lpredicted)) +  # agregamos los puntos de los valores predichos
  geom_point(color = "red") +  # agregamos los puntos de las observaciones en color rojo
  labs(title = "Grafica de residuales", x = "IMC", y = "Peso")
## `geom_smooth()` using formula 'y ~ x'

El grafico de los residuales nos permitio darnos cuanta que los valores no poseen outliers que pudieran arruinar el modelo, ademas de que la distancia de los residuos no difiere mucho entre los datos.

WCGS$index <- 1:nrow(WCGS)
ggplot(WCGS, aes(x = index, y = Lresiduals)) + 
  geom_point() +
  labs(title = "Grafica de residuales", x = "Índice", y = "Residuales")

Se logra observar bastante dispersión en los residuales lo cual da confianza respecto a nuestro modelo de regresión lineal.

Regresión múltiple

Peso

table(WCGS["chd69"])
## 
##  Ausencia Presencia 
##      2850       249
WCGS$Colesterol <- Values$Colesterol
WCGS$Edad <- Values$Edad
WCGS["chd69"][WCGS["chd69"] == "Ausencia"] <- 0
WCGS["chd69"][WCGS["chd69"] == "Presencia"] <- 1
modelo_multiple <- lm(data = WCGS , NormPeso ~ age + height + sbp + dbp + chol + behpat + bmi + dibpat + smoke )
summary(modelo_multiple)
## 
## Call:
## lm(formula = NormPeso ~ age + height + sbp + dbp + chol + behpat + 
##     bmi + dibpat + smoke, data = WCGS)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.53944 -0.03135  0.01791  0.06663  0.66852 
## 
## Coefficients: (1 not defined because of singularities)
##                Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)  -2.441e+01  7.086e-02 -344.548   <2e-16 ***
## age           3.328e-04  4.006e-04    0.831   0.4061    
## height        2.358e-01  8.605e-04  274.048   <2e-16 ***
## sbp           1.461e-04  2.266e-04    0.645   0.5190    
## dbp          -5.042e-04  3.606e-04   -1.398   0.1621    
## chol          1.138e-04  5.144e-05    2.212   0.0271 *  
## behpatA2      5.332e-03  8.138e-03    0.655   0.5124    
## behpatB3      2.762e-03  8.224e-03    0.336   0.7370    
## behpatB4      7.330e-03  9.891e-03    0.741   0.4587    
## bmi           3.240e-01  9.023e-04  359.034   <2e-16 ***
## dibpatType B         NA         NA       NA       NA    
## smokeYes     -1.032e-02  4.409e-03   -2.342   0.0193 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1194 on 3088 degrees of freedom
## Multiple R-squared:  0.9857, Adjusted R-squared:  0.9857 
## F-statistic: 2.135e+04 on 10 and 3088 DF,  p-value: < 2.2e-16
# Generar grafico multivariable
ggplot(WCGS, aes(x = bmi, y = NormPeso, color = Colesterol)) +
  geom_smooth(method = lm) +  # añadimos la linea de regresion 
  geom_point(aes(shape = as.factor(smoke))) +  # agregamos los puntos de las observaciones en color rojo
  labs(title = "Modelo lineal de variables significativas", x = "IMC", y = "Peso Normalizado", color = "Colesterol", shape = "Fuma")
## `geom_smooth()` using formula 'y ~ x'

# Generar grafico multivariable
ggplot(WCGS, aes(x = sbp, y = NormPeso, color = behpat)) +
  geom_smooth(method = lm) +  # añadimos la linea de regresion 
  geom_point(aes(shape = as.factor(Edad))) +  # agregamos los puntos de las observaciones en color rojo
  labs(title = "Modelo lineal de variables no significativas", x = "Sistolica", y = "Peso Normalizado", color = "Conducta", shape = "Edad")
## `geom_smooth()` using formula 'y ~ x'

Colesterol

Ahora que hemos buscado predecir el peso, intentaremos predecir la variable colesterol a partir de un modelo de regresión lineal multiple.

WCGS$NormChol <- bestNormalize(WCGS$chol)$x.t
head(WCGS$NormChol)
## [1]  0.5810888 -0.7282938  0.7768513 -1.2919720 -0.2271993 -0.4238293
modelo_multiple <- lm(data = WCGS , NormChol ~ age + height + sbp + dbp + behpat + bmi + dibpat + smoke )
summary(modelo_multiple)
## 
## Call:
## lm(formula = NormChol ~ age + height + sbp + dbp + behpat + bmi + 
##     dibpat + smoke, data = WCGS)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7251 -0.6251 -0.0005  0.6389  3.4185 
## 
## Coefficients: (1 not defined because of singularities)
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.5190401  0.5707442   0.909 0.363205    
## age           0.0115499  0.0032740   3.528 0.000425 ***
## height       -0.0343143  0.0070197  -4.888 1.07e-06 ***
## sbp           0.0006415  0.0018551   0.346 0.729506    
## dbp           0.0119808  0.0029448   4.068 4.85e-05 ***
## behpatA2     -0.0986521  0.0666172  -1.481 0.138741    
## behpatB3     -0.1464335  0.0672976  -2.176 0.029638 *  
## behpatB4     -0.1774567  0.0809329  -2.193 0.028408 *  
## bmi           0.0119711  0.0073856   1.621 0.105145    
## dibpatType B         NA         NA      NA       NA    
## smokeYes      0.2122971  0.0358913   5.915 3.68e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9781 on 3089 degrees of freedom
## Multiple R-squared:  0.04601,    Adjusted R-squared:  0.04323 
## F-statistic: 16.55 on 9 and 3089 DF,  p-value: < 2.2e-16
# Generar grafico multivariable
ggplot(WCGS, aes(x = dbp, y = NormChol, color = age)) +
  geom_smooth(method = lm) +  # añadimos la linea de regresion 
  geom_point(aes(shape = as.factor(smoke))) +  # agregamos los puntos de las observaciones en color rojo
  labs(title = "Modelo lineal de variables significativas", x = "DBP", y = "Colesterol Normalizado", color = "Edad", shape = "Fuma")
## `geom_smooth()` using formula 'y ~ x'

Este analisis en comparación al de la predicción de peso no parece seguir una tendencia evidente, a pesar de que el modelo lm nos indique que si hay variables significantes, el grafico nos demuestra que hacen falta mas datos y posiblemente mas variables para poder explicar los posibles niveles de colesterol de un paciente.

Regresión logística

# Verificar tipo de dato
typeof(WCGS$chd69)
## [1] "character"
head(WCGS$chd69)
## [1] "0" "0" "0" "0" "0" "0"
# Conversion numerica
WCGS$chd69 <- as.numeric(WCGS$chd69)
# Modelo logistico
modelo_logistico1 <- glm(data = WCGS, chd69 ~ bmi, family = "binomial")
summary(modelo_logistico1)
## 
## Call:
## glm(formula = chd69 ~ bmi, family = "binomial", data = WCGS)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.6778  -0.4269  -0.4004  -0.3708   2.4313  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -4.51754    0.61506  -7.345 2.06e-13 ***
## bmi          0.08402    0.02446   3.435 0.000594 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1733.1  on 3098  degrees of freedom
## Residual deviance: 1721.7  on 3097  degrees of freedom
## AIC: 1725.7
## 
## Number of Fisher Scoring iterations: 5
# Grafico de modelo logistico
ggplot(WCGS, aes(x=bmi, y=(chd69))) + 
  geom_point(alpha=.5) +
  stat_smooth(method = glm, se = FALSE, method.args = list(family = binomial)) + labs(title = "Grafica de dispersion entre Coronaria y IMC", x = "IMC", y = "Coronaria")
## `geom_smooth()` using formula 'y ~ x'

# Conversion numerica
WCGS$chd69 <- as.numeric(WCGS$chd69)
# Modelo logistico
modelo_logistico2 <- glm(data = WCGS, chd69 ~ chol, family = "binomial")
summary(modelo_logistico2)
## 
## Call:
## glm(formula = chd69 ~ chol, family = "binomial", data = WCGS)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0023  -0.4400  -0.3667  -0.3068   2.6292  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -5.320582   0.368181 -14.451   <2e-16 ***
## chol         0.012234   0.001477   8.282   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1733.1  on 3098  degrees of freedom
## Residual deviance: 1665.3  on 3097  degrees of freedom
## AIC: 1669.3
## 
## Number of Fisher Scoring iterations: 5
# Grafico de modelo logistico
ggplot(WCGS, aes(x=chol, y=(chd69))) + 
  geom_point(alpha=.5) +
  stat_smooth(method = glm, se = FALSE, method.args = list(family = binomial)) + labs(title = "Grafica de dispersion entre Coronaria y Colesterol", x = "Colesterol", y = "Coronaria")
## `geom_smooth()` using formula 'y ~ x'

Regresión logistica múltiple

modelo_logistico_mult <- glm(data = WCGS, chd69 ~  age + height + weight + sbp + dbp + chol + behpat + bmi + dibpat + smoke , family = "binomial")
summary(modelo_logistico_mult)
## 
## Call:
## glm(formula = chd69 ~ age + height + weight + sbp + dbp + chol + 
##     behpat + bmi + dibpat + smoke, family = "binomial", data = WCGS)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.2845  -0.4375  -0.3161  -0.2223   2.8600  
## 
## Coefficients: (1 not defined because of singularities)
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -24.004501  15.996764  -1.501  0.13346    
## age            0.064946   0.012295   5.282 1.28e-07 ***
## height         0.171716   0.227624   0.754  0.45062    
## weight        -0.022779   0.044936  -0.507  0.61221    
## sbp            0.018320   0.006492   2.822  0.00477 ** 
## dbp           -0.000710   0.010989  -0.065  0.94848    
## chol           0.010777   0.001569   6.870 6.41e-12 ***
## behpatA2       0.026442   0.222825   0.119  0.90554    
## behpatB3      -0.660121   0.243965  -2.706  0.00681 ** 
## behpatB4      -0.520022   0.320468  -1.623  0.10465    
## bmi            0.219331   0.314160   0.698  0.48508    
## dibpatType B         NA         NA      NA       NA    
## smokeYes       0.624085   0.143983   4.334 1.46e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1733.1  on 3098  degrees of freedom
## Residual deviance: 1549.1  on 3087  degrees of freedom
## AIC: 1573.1
## 
## Number of Fisher Scoring iterations: 6
# Grafico de modelo logistico
ggplot(WCGS, aes(x=bmi + chol + sbp + age , y=(chd69))) + 
  geom_point(alpha=.5) +
  stat_smooth(method = glm, se = FALSE, method.args = list(family = binomial)) + labs(title = "Grafica de dispersion entre Coronaria e IMC + Col + Edad + SBP", x = "IMC", y = "Coronaria")
## `geom_smooth()` using formula 'y ~ x'

No se forma la sigmoidea esperada, por lo que significa que hacen falta datos para determinar las personas que presentan una enfermedad coronaria.

0ddc52315f99eae84bc4c86d004e3e71e44e92db